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Super-cooled liquids are characterized by their fragility: the slowing 
down of the dynamics under cooling is more sudden and the jump of 
specific heat at the glass transition is generally larger in fragile liq- 
uids than in strong ones. Despite the importance of this quantity in 
classifying liquids, explaining what aspects of the microscopic struc- 
ture controls fragility remains a challenge. Surprisingly, experiments 
indicate that the linear elasticity of the glass - a purely local property 
of the free energy landscape - is a good predictor of fragility. In par- 
ticular, materials presenting a large excess of soft elastic modes, the 
so-called boson peak, are strong. This is also the case for network 
liquids near the rigidity percolation, known to affect elasticity. Here 
we introduce a model of the glass transition based on the assumption 
that particles can organize locally into distinct configurations, which 
are coupled spatially via elasticity. The model captures the men- 
tioned observations connecting elasticity and fragility. We find that 
materials presenting an abundance of soft elastic modes have little 
elastic frustration: energy is insensitive to most directions in phase 
space, leading to a small jump of specific heat. In this framework 
strong liquids turn out to lie the closest to a critical point associ- 
ated with a rigidity or jamming transition, and their thermodynamic 
properties are related to the problem of number partitioning and to 
Hopfield nets in the limit of small memory. 
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Introduction 

When a liquid is cooled rapidly to avoid crystallization, its 
viscosity increases up to the glass transition where the ma- 
terial becomes solid. Although this phenomenon was already 
used in ancient times to mold artifacts, the nature of the glass 
transition and the microscopic cause for the slowing down of 
the dynamics remain controversial. Glass-forming liquids are 
characterized by their fragility [1,2]: the least fragile liquids 
are called strong, and their characteristic time scale r fol- 
lows approximatively an Arrhenius law r(T) ^ e:>^p{Ea/kBT), 
where the activation energy Ea is independent of tempera- 
ture. Instead in fragile liquids the activation energy grows as 
the temperature decreases, leading to a sudden slowing-down 
of the dynamics. The fragility of liquids strongly correlates 
with their thermodynamic properties [3,4]: the jump in the 
specific heat that characterizes the glass transition is large in 
fragile liquids and moderate in strong ones. Various theoreti- 
cal works [5-8], starting with Adam and Gibbs, have proposed 
explanations for such correlations. By contrast few proposi- 
tions, see e.g. [9-11], have been made to understand which 
aspects of the microscopic structure of a liquid determines its 
fragility and the amplitude of the jump in the specific heat at 
the transition. 

Observations indicate that the linear elasticity of the glass 
is a key factor determining fragility - a fact a priori surpris- 
ing since linear elasticity is a local property of the energy 
landscape, whereas fragility is a non-local property charac- 
terizing transition between meta-stables states. In particu- 
lar (i) glasses are known to present an excess of soft elastic 
modes with respect to Debye vibrations at low frequencies, 
the so-called boson peak that appears in scattering measure- 
ments [12]. The amplitude of the boson peak is strongly anti- 
correlated with fragility, both in network and molecular liq- 



uids: structures presenting an abundance of soft elastic modes 
tend to be strong [13, 14]. (ii) In network glasses, where par- 
ticles interact via covalent bonds and via the much weaker 
Van der Waals interactions, the microscopic structure and 
the elasticity can be monitored by changing continuously the 
composition of compounds [15-17]. As the average valence r 
is increased toward some threshold rc, the covalent networks 
display a rigidity transition [18, 19] where the number of co- 
valent bonds is just sufficient to guarantee mechanical stabil- 
ity. Rigidity percolation has striking effects on the thermal 
properties of super-cooled liquids: in its vicinity, liquids are 
strong [20] and the jump of specific heat is small [15]; whereas 
they become fragile with a large jump in specific heat both 
when the valence is increased, and decreased [15,20]. It was 
argued [9] that fragility should decrease with valence, at least 
when the valence is small. There is no explanation however 
why increasing the valence affects the glass transition prop- 
erties in a non-monotonic way, and why such properties are 
extremal when the covalent network acquires rigidity [21]. 

Recently it has been shown that the presence of soft 
modes in various amorphous materials, including granular me- 
dia [22-25], Lennard- Jones glasses [23,26], colloidal suspen- 
sions [27-29] and silica glass [23, 30] was controlled by the 
proximity of a jamming transition [31], a sort of rigidity tran- 
sition that occurs for example when purely repulsive particles 
are decompressed toward vanishing pressure [24]. Near the 
jamming transition spatial ffuctuations play a limited role and 
simple theoretical arguments [22, 23] capture the connection 
between elasticity and structure. They imply that soft modes 
must be abundant near the transition, suggesting a link be- 
tween observations (i) and (ii). However these results apply 
to linear elasticity and cannot explain intrinsically non-linear 
phenomena such as those governing fragility or the jump of 
specific heat. In this article we propose to bridge that gap 
by introducing a model for the structural relaxation in super- 
cooled liquids. Our starting assumption is that particles can 
organize locally into distinct configurations, which are cou- 
pled at different points in space via elasticity. We study what 
is perhaps the simplest model realizing this idea, and show 
numerically that it captures qualitatively the relationships be- 
tween elasticity, rigidity, thermodynamics and fragility. The 
thermodynamic properties of this model can be treated theo- 
retically within a good accuracy in the temperature range we 
explore. Our key result is the following physical picture: when 
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Fig. 1. Top row: sketches of covalent networks with different mean valence r 
around the valence Vc'- red solid lines represent covalent bonds; cyan dash lines rep- 
resent van der Waals interactions. Bottom: sketch of our elastic network model with 
varying coordination number z (defined as the average number of strong springs in 
red) around Maxwell threshold Zc\ cyan springs have a much weaker stiffness, and 
model weak interactions. 

there is an abundance of soft elastic modes, elastic frustration 
vanishes, in the sense that a limited number of directions in 
phase space cost energy. Only those directions contribute to 
the specific heat, which is thus small. Away from the critical 
point, elastic frustration increases: more degrees of freedom 
contribute to the jump of specific heat, which increases while 
the boson peak is reduced. 

Model 

Our main assumption is that in a super-cooled liquid, nearby 
particles can organize themselves into a few distinct configu- 
rations. Consider for example covalent networks sketched in 
Fig. 1, where we use the label {ij) to indicate the existence 
of a covalent bond between particles i and j. If two covalent 
bonds (12) and (34) are adjacent, there exists locally another 
configuration for which these bonds are broken, and where the 
bonds (13) and (24) are formed instead. These two configura- 
tions do not have the same energy in general. Moreover going 
from one configuration to the other generates a local strain, 
which creates an elastic stress that propagates in space. In 
turn, this stress changes the energy difference between local 
configurations elsewhere in the system. This process leads to 
an effective interaction between local configurations at differ- 
ent locations. 

Our contention is that even a simple description of the 
local configurations - in our case we will consider two-level 
systems, and we will make the approximation that the elastic 
properties do not depend on the levels - can capture several 
unexplained aspects of super-cooled liquids, as long as the 
salient features of the elasticity of amorphous materials are 
taken into account. To incorporate in particular the pres- 
ence of soft modes in the vibrational spectrum we consider 
random elastic networks. The elasticity of three types of net- 
works have been studied extensively: networks of springs ran- 
domly deposited on a lattice [32] , on- lattice self-organized net- 
works [33] and off-lattice random networks with small spatial 
fluctuations of coordination [34-36]. We shall consider the 
last class of networks, which are known to capture correctly 



the scaling properties of elasticity near jamming, and can be 
treated analytically [22,35,37]. It is straightforward to ex- 
tend our model to self-organized networks ^. In our model 
two kinds of springs connect the N nodes of the network: 
strong ones, of stiffness k and coordination and weak ones, 
of stiffness and coordination z^. These networks undergo 
a rigidity transition as z crosses Zc = 2d, where d is the spa- 
tial dimension. For z < Zc elastic stability is guaranteed by 
the presence of the weak springs. As indicated in Fig. 1, this 
situation is similar to covalent networks, where the weak Van 
der Waals interactions are required to insure stability when 
the valence r is smaller than its critical value Vc- 

Initially when our network is built, every spring {ij) is at 
rest: the rest length follows = | |R° — R° 1 1 , where R° is the 
initial position of the node i. To allow for local changes of con- 
figurations we shall consider that any strong spring {ij) can 
switch between two rest lengths: = + ea(^ij), where 

cr^^ij) = ±1 is a spin variable. There are thus two types of vari- 
ables: the Ns = zN/2 spin variables {o-{ij)}^ which we shall 
denote using the ket notation |cr), and the Nd coordinates of 
the particles denoted by |R). The elastic energy ^(|R), \cr)) is 
a function of both types of variables. The inherent structure 
energy 7-i{\a)) associated with any configuration \a) is defined 
as: 

n{\a)) = min|R)^(|R), \a)) = ke^H{\cj)) [1] 

where we have introduced the dimensionless Hamiltonian 
We shall consider the limit of small e, where the vibrational en- 
ergy is simply that of harmonic oscillators. In this limit all the 
relevant information is contained in the inherent structures 
energy, since including the vibrational energy would increase 
the specific heat by a constant, which does not contribute to 
the jump of that quantity at the glass transition. In this limit, 
linear elasticity implies the form: 

n{\a)) = \Y.Q,,pcj,ap + o(e^) ^ \{cj\g\cj) + o{^) [2] 

where 7 and (3 label strong springs, Q-f,p is the Green func- 
tion describing how a dipole of force applied on the contact 7 
changes the amplitude of the force in the contact (3. Note that 
models where some kind of defects interact elastically, lead- 
ing to Hamiltonians similar in spirit to that of Eq.(2), have 
been proposed to investigate the low-temperature properties 
of glasses [38] and super-cooled liquids [10]. These models 
however assume continuum elasticity, unlike our model which 
can incorporate the effects of a rigidity transition and the 
presence of a boson peak. 

^^,/3 is computed in Supplementary Information (SI) part 
A and reads: 

g ^X-SsM-^Sl [3] 

where T is the identity matrix, and Ss and the dimensionless 
stiffness matrix M are standard linear operators connecting 
forces and displacements in elastic networks [39]. They can 
be formally written as: 

M = SlS. + ^SiS^, [4] 
k 

s, = 

where (i|R) = R^, (u)» indicate a summation over the strong 
springs (• = s) or the weak springs (• = w). <S, is a A/', x dN 



^If the network is assumed to present a rigidity window for which a subpart 
of the system is critical, we expect that within this window thermodynamics 
and fragility would behave as the critical point 2 = 2^ in the present model. 
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matrix which projects any displacement field onto the contact 
space of strong or weak springs. The components of this lin- 
ear operator are uniquely determined by the unit vectors Uij 
directed along the contacts (ij) and point toward the node i. 

Finally note that the topology of the elastic network is 
frozen in our model. This addition of frozen disorder is obvi- 
ously an approximation, as the topology itself should evolve as 
local configurations change. Our model thus misses the evo- 
lution of elasticity with temperature, that presumably affects 
the slowing down of the dynamics [8] and gives a vibrational 
contribution to the jump of specific heat [8,40]. Building 
models which incorporate this possibility, while still tractable 
numerically and theoretically, remains a challenge. 
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Fig. 2. Specific heat c(T) v.s. rescaled temperature T/Tg for various excess 
coordination Sz = z — Zc indicated in legend, for o; = 3 X 10~^ and d = 2. 
c(T) displays a jump at the glass transition. Solid lines are theoretical predictions, 
deprived of any fitting parameters, of our mean-field approximation. They terminate 
at the Kautzman temperature T^. Inset: glass transition temperature Tg vs Sz for 
several amplitude of weak interactions a, as indicated in legend. 
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Fig. 3. Jump of specific heat Ac versus excess coordination Sz \n d = 2 for 
different strength of weak springs a, as indicated in legend. Solid lines are mean-field 
predictions not enforcing the orthogonality of the \Srp), dashed-line corresponds to 
the ROM where orthogonality is enforced. In both cases the specific heat is computed 
at the numerically obtained temperature Tg. Inset: theoretical predictions for Ac 
vs Sz computed at the theoretical temperature T^. 




Simulation 

Network structure. Random networks with weak spatial fluc- 
tuations of coordination can be generated from random pack- 
ings of compressed soft particles [34-36]. We consider pack- 
ings with periodic boundary conditions. The centers of the 
particles correspond to the nodes of the network, of unit mass 
m = 1, and un-stretched springs of stiffness /c = 1 are put be- 
tween particles in contact. Then springs are removed, prefer- 
ably where the local coordination is high, so as to achieve the 
desired coordination z. In a second phase, A^w weak springs 
are added between the closest unconnected pairs of nodes. 
The relative effect of those weak springs is best characterized 
hy a = {z^/d){k^/k), which we modulate by fixing = 6 
and changing {k.^/k). Note that an order of magnitude esti- 
mate of a in covalent glasses can be obtained by comparing the 
behavior of the shear modulus G in the elastic networks [34] 
and in network glasses near the rigidity transition. As shown 
In Fig. 6 of SI part B, this comparison yields the estimate 
that a e [0.01,0.05]. 

Thermodynamics. We introduce the rescaled temperature T = 
TkB / (ke^) where T is the temperature. To equilibrate the 
system, we perform a one spin-flip Monte Carlo algorithm. 
The energy Ti of configurations are computed using Eq.(2). 
We use 5 networks oi N — 256 nodes in two dimensions and 
N — 216 in three dimensions, each run with 10 different initial 
configurations. Thus our results are averaged on these 50 re- 
alizations. We perform 10^ Monte Carlo steps at each T. The 
time-average inherent structure energy E = {V) is calculated 
as a function of temperature, together with the specific heat 
Cv = dE/dT. The intensive quantity c(T) = Cv/Ns is repre- 
sented in Fig. 2 for various excess coordination 6z = z — Zc 
and a = 3 X 10~^. We observe that the specific heat increases 
under cooling, until the glass transition temperature Tg where 
c(T) rapidly vanishes, indicating that the system falls out of 
equilibrium. 

The amplitude of c(T) just above Tg thus corresponds 
to the jump of specific heat Ac ^, and is shown in Fig. 3. 
Our key finding is that as the coordination increases, Ac{z) 
varies non-monotonically and is minimal in the vicinity of the 
rigidity transition for all values of a investigated, as observed 
experimentally [15]. This behavior appears to result from a 
sharp asymmetric transition at a ^ 0. For z > Zc we observe 
that Ac{z) oc 6z. The jump in specific heat thus vanishes as 
6z ^ where the system can be called "perfectly strong". 
For z < Zc, Ac is very rapidly of order one. When a increases, 
this sharp transition becomes a cross-over, marked by a min- 
imum of Ac(^) at some coordination larger but close to Zc. 

Dynamics. To characterize the dynamics we compute the cor- 
relation function C{t) = ^{a{t)\a{0)), which decays to zero 
at long time in the liquid phase. We define the relaxation 
time r as C(r) = 1/2, and the glass transition temperature 
Tg as r{Tg) /t[oo) — 10^. Finite size effects on r appear to 
be weak, as shown in SI part C. The Angell plot represent- 
ing the dependence of r with inverse rescaled temperature is 
shown in the inset of Fig. 4. It is found that the dynamics 



In our approach the absolute value of the specific heat will depend on the 
minimal number of particles needed to generate distinct local configurations, 
and on the number of configurations such a group can generate. If those two 
numbers are of order one, our model predicts that Ac is of order one per 
particle or "beads", as observed near the glass transition [6]. 
^When q; — )■ and 5z > 0, (^bp ^^'i D{ujbp) 1 [22], whereas 

G ^ Sz [23], leading to 1/Abp ~ VSz- For 5z < 0, G ^ -a/5z [34]. On the 
other hand the boson peak is governed by the fraction ~ Sz of floppy modes, 
which gain a finite frequency ~ ^/o^ [35] thus we expect D^ujj^p) ~ —5z/-s/ol 
and 1/Abp ~ y/—5z. 
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follows an Arrhenius behavior for a ^ and z ^ Zc. Away 
from the rigidity transition, the slowing down of the dynamics 
is faster than Arrhenius. To quantify this effect we compute 
the fragility m = ^j^^ff^\T=Tg , whose variation with coordi- 
nation is presented in Fig. 4. Our key finding is that for all 
weak interaction amplitudes a studied, the fragility depends 
non-monotonically on coordination and is minimal near the 
rigidity transition, again as observed empirically [20] in cova- 
lent liquids. As was the case for the thermodynamic proper- 
ties, the fragility appears to be controlled by a critical point 
present at a = and z — Zc where the liquid is strong, and the 
dynamics is simply Arrhenius. As the coordination changes 
and \bz\ increases, the liquid becomes more fragile. The rapid 
change of fragility near the rigidity transition is smoothed over 
when the amplitude of the weak interaction a is increased. 



Correlating boson peak and fragility. The presence of soft elas- 
tic modes in glasses can be analyzed by considering the max- 
imum of Z(uS) = B^ilj) \y^^ where D(uS) is the vibrational 
density of states and Dd{(^) oc cj^ the Debye model for this 
quantity. The maximum of Z{oj), Abp = Z{ujbp) defines the 
boson peak frequency ujbp [12] and it normalized amplitude 
Abp [13]. The inverse oi Abp was shown to strongly correlate 
with fragility [13, 14] both in molecular liquids and covalent 
networks. 

To test if our model can capture this behavior we com- 
pute the density of states via a direct diagonalization of the 
stiffness matrix, see Eq.(4). To compute Z{ijj) the Debye den- 
sity of states is estimated as D(ijo) ^ lkP' jG^^^ where G is 
the shear modulus (bulk and shear moduli scale identically 
in this model, see e.g. [37]). Then we extract the maximum 
Abp — Z{(jJbp)- The dependence of 1/Abp is represented 
in Fig. 5 and shows a minimum near the rigidity transition, 
and even a cusp in the limit a ^ 0. This behavior can be 
explained in terms of previous theoretical results on the den- 
sity of states near the rigidity transition, that supports that 
1/Abp - \Sz\^^^ when a^O^. 

Fig. 5 shows that 1/Abp and the liquid fragility m are cor- 
related in our model, thus capturing observations in molecular 
liquids. The model also predicts that 1/Abp and the jump 
of specific heat are correlated. Note that the correlation be- 
tween fragility and Abp is not perfect, and that two branches, 
for glasses with low and with high coordinations, are clearly 
distinguishable. In general we expect physical properties to 
depend on the full structure of the density of states, as will 
be made clear for the thermodynamics of our model below. 
The variable Abp, which is a single number, cannot capture 
fully this relationship. In our framework it is a useful quantity 
however, as it characterizes well the proximity of the jamming 
transition. 



Theory 

Thermodynamics in the absence of weak interactions {(X = 

0). In the absence of weak springs the thermodynamics is non- 
trivial il z > Zc, otherwise the inherent structure energies are 
all zero. Then Eq.(4) implies M = S^Ss, and Eq.(3) leads to 
Q — T — Ss{SlSs)~^ Sl . Inspection of this expression indicates 
that ^ is a projector on the kernel of <Ss, which is generically of 
dimension Ns — Nd = 6zN/2. This kernel corresponds to all 
the sets of contact forces that balance forces on each node [23] . 
We denote by \Srp),p = 1, ...,SzN/2 an orthonormal basis of 
this space. We may then rewrite Q = Eq.(2) 



SzN/2 



[5] 



Eq.(5) is a key result, as it implies that near the rigidity tran- 
sition the number SzN/2 of directions of phase space that cost 
energy vanishes. Only those directions can contribute to the 
specific heat, which must thus vanish linearly in 6z as the 
rigidity transition is approached from above. 

Eq.(5) also makes a connection between strong liquids in 
our framework and well-know problems in statistical mechan- 
ics. In particular Eq.(5) is similar to that describing Hop- 
field nets [41] used to store SzN/2 memories consisting of the 
spin states \Srp). The key difference is the sign: in Hopfield 
nets memories correspond to met a- stables states, whereas in 
our model the vectors \Srp) corresponds to maxima of the 
energy. A particularly interesting case is 6zN/2 = 1, the 
closest point to the jamming transition which is non-trivial. 
In this situation the sum in Eq.(5) contains only one term: 

n(\a)) = ^((7\6rif = |(E^=i^n,c.cra)^- This Hamiltonian 
corresponds to the NP complete partitioning problem [42], 
where given a list of numbers (the Sri^a) one must partition 
this list into two groups whose sums are as identical as pos- 
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Fig. 4. Fragility m and rescaled fragility nisc extrapolated to the experimental 
dynamical range (see SI part D for a definition) versus excess coordination 6z for 
different strength of weak interactions a as indicated in legend, In d = 2. Dash 
dot lines are guide to the eyes, and reveal the non-monotonic behavior of m near the 
rigidity transition. Inset: Angell plot representing logr v.s. inverse temperature 
Tg/T for different 6z and a = 3 X 10~^. 
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Fig. 5. Left: I nverse boson peak amplitude 1/Abp versus excess coordination 
6z in our d = 3 elastic network model, for different weak interaction strenghts as 
indicated in legend. Dash dot lines are drawn to guide one's eyes. Right: Inverse 
boson peak amplitude 1/Abp versus fragility m for different weak springs a. 
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sible. Thermodynamically this problem is known [43] to map 
into the random energy model [44] where energy levels are 
randomly distributed. 

It is in general very difficult to compute the thermody- 
namic functions of the problem defined by Eq.(5) because the 
vectors \5rp) present spatial correlations, as must be the case 
since the amplitude of the interaction kernel Q^^p must decay 
with distance. However this effect is expected to be mild near 
the rigidity transition. Indeed there exists a diverging length 
scale at the transition, see [35] for a recent discussion, below 
which Q^^p is dominated by fluctuations and decays mildly 
with distance. Beyond this length scale Q^^p presents a dipo- 
lar structure, as in a standard continuous elastic medium. We 
shall thus assume that \Srp) are random unitary vectors, an 
approximation of mean-field character expected to be good 
near the rigidity transition. 

Within this approximation, the thermodynamic properties 
can be derived for any spectrum of Q [45] . If the orthogonal- 
ity of the vectors \5rp) is preserved, the Hamiltonian of Eq.(5) 
corresponds to the Random Orthogonal Model (ROM) whose 
thermodynamic properties have been derived [45] as well as 
some aspects of the dynamics [46] . Comparison of the specific 
heat of our model and the ROM predictions of [45] is shown in 
Fig. 3 and are found to be very similar. For sake of simplicity, 
in what follows we shall also relax the orthogonal condition 
on the vectors \Srp). This approximation allows for a straight- 
forward analytical treatment in the general case a 7^ 0, and is 
also very accurate near the rigidity transition since the number 
of vectors SzN/2 is significantly smaller than the dimension of 
the space dN they live in, making random vectors effectively 
orthogonal. Under these assumptions we recover the random 
Hopfield model with negative temperature. 

In the parameter range of interest, the Hopfield free en- 
ergy — ln(.2^) (here (...) represents the disorder average 
on the \Srp)) is approximated very precisely by the annealed 
free energy ln(Z) (this is obviously true for the number par- 
titioning problem that maps to the Random Energy Model), 
which can be easily calculated. Indeed in our approximations 
the quantities Xy = {a\Srp) are independent gaussian random 
variables of variance one, and: 



/ NSz/2 



Z oc 



/ n 

^ \ p=l 



27r 



Performing the Gaussian integrals we find: 
c(T) ' 



2z (1 + T)2 



[6] 



[7] 



The Kautzman temperature defined as s{Tk) = is found 
to follow Tk ^ |2"^'=/^'=. Eq.(7) evaluated at Tg is tested 
against the numerics in Fig. 3 and performs remarkably well 
for the range of coordination probed. 

General case (a 7^ 0).To solve our model analytically in 
the presence of weak interactions, we make the additional ap- 
proximation that the associated coordination ^ 00, while 
keeping a = z^k^/{kd) constant. In this limit weak springs 
lead to an effective interaction between each node and the 
center of mass of the system, that is motionless. Thus the 
restoring force stemming from weak interactions |Fw) follows 
|Fw) = — a|^R), leading to a simple expression in the stiffness 
matrix Eq.(4) for the weak spring contribution ^Sl^Sw = olL. 
It is useful to perform the eigenvalue decomposition: 



where |(5Ra;) is the vibrational mode of frequency uj in the 
elastic network without weak interactions. We introduce the 
orthonormal eigenvectors in contact space l^r^;) = «Ss|^Ru;)/^ 
defined for > 0. For < these vectors form a complete 
basis of that space, of dimension Ns- When 8z > ^ however, 
this set is of dimension Nd < Ns, and it must be completed 
by the kernel of tSg, i.e. the set of the \5rp),p = 1, ...,5zN/2 
previously introduced. Using this decomposition in Eq.(3,4) 
we find: 

SzN/2 

p=l oj>0 

where the first term exists only for 6z > 0. Using the mean 
field approximation that the set of \Srp) and \Sruj) are random 
gaussian vectors, the annealed free energy is readily computed, 
as shown in SI part E. We find in particular for the specific 
heat: 



c{T) = 



6z 0{8z) 1 
2^(1 + T)2 ^ 2N, 



a + (w2 ^ a)T 



fiol 



where 0{x) is the unitary step function. To compare this pre- 
diction with our numerics without fitting parameters, we com- 
pute numerically the vibrational frequencies for each value of 
the coordination. Our results are again in excellent agreement 
with our observations, as appears in Figs. 5, 5. 

To obtain the asymptotic behavior near jamming, we re- 
place the summation over frequencies in Eq.(lO) by an integral 
Sa;>o ~^ / d^D(ct;). The associated density of vibrational 
modes D((jo) in such networks has been computed theoreti- 
cally [22,35,37]. These results allows us to compute the scal- 
ing behavior of thermodynamic properties near the rigidity 
transition, see SI part F. We find that the specific heat in- 
creases monotonically with decreasing temperature. Its value 
at the Kautzman temperature thus yields an upperbound on 
the jump of specific heat. In the limit a ^ 0, we find that 
a sudden discontinuity of the jump of specific heat occurs at 
the rigidity transition: 



6z 
2z 



lim c{Tk) ~ — 

5z^0- Sz 



for 6z > 



[11] 

\12] 



[8] 



Eq.(ll) states that adding weak interactions is not a singular 
perturbation for 5z > 0, and we recover Eq.(7). On the other 
hand for Sz < 0, the energy of inherent structures is zero in 
the absence of weak springs, which thus have a singular ef- 
fect. The relevant scale of temperature is then a function of 
a. In particular we find that the Kautzman temperature is 
sufficiently low that all the terms in the second sum of Eq.(lO) 
contribute significantly to the specific heat, which is therefore 
large as Eq.(12) implies. Thus as the coordination decreases 
below the rigidity transition, one goes discontinuously from 
a regime where at the relevant temperature scale the energy 
landscape consists of a vanishing number of costly directions 
in phase space, whose cost is governed by the strong interac- 
tion /c, to a regime where the weak interaction a is the relevant 
one, and where at the relevant temperature scale all directions 
in phase space contribute to the specific heat. 

Note that although the sharp change of thermodynamic 
behavior that occurs at the rigidity transition is important 
conceptually, empirically a smooth cross-over will always be 
observed. This is the case because (i) a is small but finite. 
As a increases this sharp discontinuity is replaced by a cross- 
over at a coordination 6z ^ ln(l/a)~^ (see SI part F) where 
c{Tk,z) is minimal, as indicated in the inset of Fig. 3. (ii) 
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The Kautzman temperature range is not accessible dynami- 
cally, i.e. Tg » Tk near the rigidity transition. Comparing 
Fig. 3 with its inset, our theory predicts that the minimum 
of c{Tg) is closer to Zc and more pronounced than at Tk- 



Discussion 

Previous work [31] has shown that well-coordinated glasses 
must have a small boson peak, which increases as the coordi- 
nation (or valence for network glasses) is decreased toward the 
jamming (or rigidity) transition. Here we have argued that as 
this process occurs, elastic frustration vanishes: thanks to the 
abundance of soft modes, any configuration (conceived here 
as a set of local arrangements of the particles) can relax more 
and more of its energy as jamming is approached from above. 
As a result, the effective number of degrees of freedom that 
cost energy and contribute to the jump of specific heat at the 
glass transition vanishes. As the coordination is decreased 
further below the rigidity transition, the scale of energy be- 
comes governed by the weak interactions (such as Van der 
Waals) responsible for the finite elasticity of the glass. At 
that scale, all direction in phase space have a significant cost 
and the specific heat increases. This view potentially explains 
why linear elasticity strongly correlates to key aspects of the 
energy landscape in network and molecular glasses [13-16]. 
This connection we propose between structure and dynam- 
ics can also be tested numerically. Our results are consistent 
with the observation that soft particles become more fragile 
when compressed away from the random close packing [47] ^. 
Another interesting parameter to manipulate is the amplitude 
of weak interactions, which can be increased by adding long- 
range forces to the interaction potential [23,26]. According to 
our analysis, doing so should increase fragility, in agreement 
with existing observations [48]. 

The model of the glass transition we introduced turns out 
to be a spin glass model, with the specificity that (i) the in- 
teraction is dipolar in the far field, and that (ii) the sign of 
the interaction is approximatively random below some length 
scale Ic that diverges near jamming, where the coupling matrix 
has a vanishingly small rank. Applying spin glass models to 
structural glasses have a long history. In particular the Ran- 
dom First Order Transition theory (RFOT) [6] is based on 
mean-field spin glass models that display a thermodynamic 
transition at some Tk where the entropy vanishes. A phe- 
nomenological description of relaxation in liquids near Tk 
based on the nucleation of random configurations leads to 
a diverging time scale and length scale ^ at Tk [6,7]. One 
limitation of this approach is that no finite dimensional spin 
models with two-body interactions have been shown to follow 
this scenario so far [49], and it would thus be important to 
know if our model does display a critical point at finite tem- 
perature. Our model will also allow one to investigate the 
generally neglected role of the action at a distance allowed by 
elasticity, characterized by a scale Ic. In super-cooled liquids 
heterogeneities of elasticity (that correlates to irreversible re- 
arrangements) can be rather extended [50] suggesting that Ic 
is large. This length scale may thus play an important role in 
a description of relaxation in liquids, and in deciphering the 
relationship between elastic and dynamical heterogeneities. 
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SUPPLEMENTARY INFORMATION 
A. Stiffness and Coupling matrices 

Consider a network of N nodes connected by Nc springs. If 
an infinitesimal displacement field |(5R) is imposed on the 
nodes, the change of length of the springs can be written 
as a vector \5r) of dimension Nc. For small displacements 
this relation is approximately linear: \5r) — <S|^R), where S 
is an Nc x Nd matrix. To simplify the notation, we write 
tS as an A^c X A/" matrix of components of dimensions 
which gives S^^i = dr^/dHi = ^^,in^, where 6j,i is non- 
zero only if the contact 7 includes the particle i, and is 
the unit vector in the direction of the contact 7, pointing to- 
ward the node i. Using the bra-ket notation, we can rewrite 
S = J2{ij)=jh)^7i{'^\ ~ 01)5 where the sum is over all the 
springs of the network. Note that the transpose <S* of S re- 
lates the set of contact forces |/) to the set |F) of unbalanced 
forces on the nodes: |F) = S^lf), which simply follows from 
the fact that = ^i,if-f''^i = J2j fi^i,i [39]. 

The stiffness matrix is a linear operator connecting 
external forces to the displacements: A4|^R) = |F). Intro- 
ducing the Nc X Nc diagonal matrix /C, whose components are 
the spring stiffnesses K,^^ = /c^, we have for harmonic springs 
I/) = lC\5r). Applying on each side of this equation, we 
get |F) = <S*|/) = <S*/C<S|^R), which thus implies [39]: 

M = S'KS. 

Let us assume that starting from a configuration where all 
springs are at rest, the rest lengths of the springs are changed 
by some amount \y). This will generate an unbalanced force 
field |F) = S^JC\y) on the nodes, leading to a displacement 
j^R) = M-^S'JCly). The elastic energy ^ = ^{y-6r\K:\y-6r) 
is minimal for this displacement and the corresponding energy 
n is: 

n\y)) = l(y\^- JCSM-'s'Jc\y). [si-i] 




8z/z^ 

Fig. 6. Squares show the shear modulus G normalized by its value at the rigidity 
threshold for Ge-Se, taken from Ref. [53]. Circles show G for Ge-Sb-Se, taken from 
Ref. [54]. Lines display the shear modulus G for network models in fi = 3 using 
different a, as indicated in the legend. 



A quantitative treatment of the vicinity of random close packing would re- 
quire including the dependence of elasticity on temperature neglected here. 
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In our model, 2/7 = for weak springs and — ea^ for strong 
springs of stiffness k, implying that JC\y) — k\y). Introducing 
the dimensionless stiffness matrix M = M/k and the restric- 
tion Sl of the operator <S* on the subspace of strong contacts 
of dimension Ns^ i.e. SI\cf) = S^\y), Eq.(SI-l) yields: 



'HiW)) = -{a\g\a) where g 



where X is the identity matrix, and Q is the coupling matrix 
used in the main text. Note that in our model the diagonal 
matrix JC contains only two types of coefficients k^ and k, cor- 
responding to the stiffnesses of weak springs and stiff springs 
respectively. Then the dimensionless stiffness matrix can be 
written as = S^Ss + ^ <S4*Sw, where <S4 is the projection 
of the operator on the subspace of weak contacts. 



B. Shear Modulus of the random elastic networks 

An explicit expression for the shear modulus G of an elastic 
network can be found using linear response theory [51,52]. In 
particular, let us consider the shear on the (x, y) plane. In the 
contact vector space, a shear strain can be written as S\y^^) 
with ^ ^ 1 which represents the amplitude of the strain and 
\y^^) corresponds to a unit shear strain. The components of 
\y^^) are given by y^^ — , where is the rest length of 

the spring 7, and Ax^ and /\y^ are its projections along the 
X and y directions. From the last section (A), we can obtain 
the total energy induced by a shear strain i-L(\y^^))5'^] hence 
the shear modulus G = 2n{\y'^))/V. 

To estimate the value of a, we consider the dependence of 
the shear modulus with coordination or valence in the vicin- 
ity of the rigidity transition, which is smooth for large a and 
sudden for small a in our networks, see Fig. 6. Compar- 
ing networks and real chalcogenide glasses we find that the 
cross-over in the elastic modulus is qualitatively reproduced 
for a e [0.01,0.05]. 



C. Finite size effects on fragility 

To estimate the role of finite size effects on the dynamics, we 
use two different system sizes N — 64 and N — 256. As 
shown in Fig. 7, the Angell plot for the relaxation time, and 
therefore our estimation of the fragility, appears to be nearly 
independent of the system size. Note, however, that the cor- 
relation function C (t) shows some finite size effects very close 
to the isostatic point = Zc, = 0), but that it does not af- 
fect our measure of r significantly. In particular we find that 
near isostaticity, the distribution of relaxation time is broad 
for small systems, and becomes less and less so when the sys- 
tem size increases. We noticed that this effect also disappears 
if a two-spin flips Monte-Carlo is used, instead of the one-spin 
flip algorithm we perform. 



D. Rescale fragility with different dynamical range 

The value of fragility depends on the definition of glass tran- 
sition, in particular on the dynamical range. In super-cooled 
liquids the glass transition occurs when the relaxation time 
is about 10^^ larger than the relaxation time at high tem- 
perature. Thus the dynamical range in experiments (which 
corresponds to the fragility of a perfectly Arrhenius liquid) is 
= 16. In our simulation, the same quantity is R — 5. It is 
possible however to rescale our values of fragility to compare 
with experimental data, if we extrapolate the dynamics. We 
shall assume a Vogel-Fulcher-Tammann (VFT) relation at low 
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Fig. 7. Angell plot representing logr v.s. inverse temperature Tg/T for differ- 
ent 5z and two system sizes = 64 and = 256, a = 0.0003. 



temperature. 



logi 



t(T) 



A 



TO T 

We define the dynamical range as: 



R = logi 



r{Tg 

To 



A 



Thus we can express the fragility as: 
91ogior(T)/ro 



niR 



dTJT 



R + R 



A • 



[SI-2] 



To and A are assumed to be independent of dynamical range. 
Using the notation m — ms and rrisc — miQ we get from 
Eq.(SI-2): 



16^ 

16+ ^(m- 



5) = 10.24m - 35.2 



The amplitude of fragility we find turn out to be comparable 
to experiments when a — 0.03, in particular for z > Zc. For 
the smallest coordination explored our results underestimate 
somewhat the fragility, slightly above 50 in our model and 
about 80 experimentally. This is not surprising considering 
that our model is phenomenological, and the extrapolation 
we made to compare different dynamical ranges. 



E. Canonical ensemble approach with weak springs 

In the case where a 7^ 0, the annealed free energy can be eas- 
ily calculated under the assumption that \5rp) and \5r^) are 
random Gaussian vectors. The Hamiltonian in Eq.(lO) can be 
rewritten as: 



1 



E 



^l...SzN/2 



where Xp — {5rp\a) and Xoj — {5roj\(T) represent independent 
random variables for each configuration |cr). In the thermo- 
dynamic limit the random variables Xp and Xoj are Gaussian 
distributed with zero mean and unit variance. The averaged 
partition function is given by: 



Footline Author 



PNAS I Issue Date | Volume | Issue Number | 7 



Z = 2" 



e 2 



27r 



^^^1^ / X -1/2 / 

=2- n 



ci;>0 



-1/2 



It turns out that the Debye regime contribution to the in- 
tegrals is neghgible near the jamming threshold. To capture 
the scaling behavior near jamming, we approximate by 
a square function. This simplified description allows further 
analytical progress while preserving the same qualitative be- 
havior. Since the Debye regime can be neglected, we choose: 



From the average partition function the density of free en- 
ergy per spring /(T) and any other thermodynamic quantities 
are readily computed. In particular, the energy density ^(T), 
the specific heat c(T) and the entropy density s(T) write: 



2Ns 



SzN/2 



e(T) 



2Ns 



a;>0 



a/T 



-Tin 2 



5zN/2 



^ 1 +T ^ ^ a 



aT 



uj>0 



c{T) 



1 

27V: 



5zN/2 ^ fa 

p=i v-^ ^ ; cu>o ^ ^ 



[SI-2] 



5zN/2 



ln(l + 



p=i 

lnfl + 



a/T 



l + T 



In the limit a ^ 0, for any finite temperature T, the sum 
over the vibration modes {uo > 0) vanishes, and we recover 
the expressions in the absence of weak springs for pure rigid 
networks. Note that Eq.(SI-2) corresponds to Eq.(ll) in the 
article. 



F. Continuous density of states limit: Analytical results 

In the thermodynamic limit N ^ oo, we can replace the sum 
over frequencies by an integral: X]^>o ~^ I dujD(u) for 
Sz < 0, and T^^^q ^dj dujD{uj) for 5z > 0. The den- 
sity of states D{uj) is the distribution of vibrational modes of 
random elastic networks, which has been computed theoreti- 
cally [22,35,37]. There are two frequency scales in the random 
network : cj* ^ \6z\ above which a plateau of soft modes ex- 
ist, and a cut-off frequency ojc I- Below cj*, rigid networks 
show plane wave modes [22, 34, 37] with a characteristic De- 
bye regime D(uj) ~ cj^, unlike fioppy networks, which show 
no modes in this gap [35]. 



Considering cj* = ^^c, the cut-off frequency cjc 1 is 
the only fitting parameter of the simplified continuum model. 
Rescaling as a — aUc, we obtain that all the thermodynamic 
functions depend uniquely on o; = ^ 5z. In partic- 

ular, the specific heat is: 



UJ* < UJ < UJc 5z <{) 

0<uj <ujc 6z> 0. 



c{T, 6z, a) = < 



^a(l + l/T) 
4z (1 + T)2 



arctan( 



7a(l + l/T) 



+ i+a(i+i/T) + arctan( 



5z/zc 



+ 



Va(l + 1/T) 
zc Sz-^/z-^+cxil + l/T) 

arctan( 



V«(l + 1/T) 



Va(l + 1/T) 
4:z (1+T)2 



+ 



Va(l + 1/T) 
l+a(l + l/T) 



2z (1+T)2 



Sz<0 



Sz > 0. 



We compute the jump of specific heat at the Kautzman 
temperature, where the entropy vanishes s{TK,Sz,a) = 0. In 
the continuous limit, the equations for Tk can be approxi- 
mated by: 



ln2 ; 



2z 



ln(l + ^) + 2^ 



- arctan i 



6z 
'2z 



0{5z) IuTk, 



where the conditions 1^2; | ^ Zc and a ^ k^/k <^ 1 have been 
used. There is no simple analytical expression for Tk, how- 
ever, one can observe the existence of two asymptotic regimes: 
Tk r^aioi 6z < |l/lnQ;| and 2"^^/^^ for 6z > |l/lna|. 

Then the specific heat at the transition temperature is given 
by: 



c(Tk, Sz, a) ~ 



5z 
2z 



6z < |l/lnQ;| 
Sz > |l/lna| 



From these asymptotic behaviors one gets that the specific 
heat display a non-monotonous behavior with coordination, 
with a minimum whose position scales as Sz ~ |l/lna|. 
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